!----------------------------------------------------------------------!
! Optional advection scheme for CCATT-BRAMS/BRAMS models version 4.2+  ! 
! Based on Walcek, 2000 (JGR) and Walcek and Aleksic, 1998 (ATENV).    ! 
! The scheme is highly conservative, monotonic and keeps mass mixing   ! 
! ratio positive definite. 					       ! 
! Implemented by Saulo Freitas (saulo.freitas@cptec.inpe.br) @ Jun/2009! 
!----------------------------------------------------------------------!

MODULE monotonic_adv
  
  USE node_mod, only:        &
          	      ibcon, &  !INTENT(IN)
          	      mynum, &  !INTENT(IN)
          	      nodei0,   &  !INTENT(IN)
          	      nodej0,   &  !INTENT(IN)
          	      nodemyp,  &  !INTENT(IN)
          	      nodemxp,  &  !INTENT(IN)
          	      nodemzp   	!INTENT(IN)
  USE mem_grid, ONLY:        &
  	    	     dtlt,   & !INTENT(IN)
  	    	     ngrids, & !INTENT(IN)
  	    	     ngrid,  & !INTENT(IN)
	    	     dzt,    & !INTENT(IN)
	    	     dztn,   & !INTENT(IN)
            	     grid_g, & !INTENT(IN)	 
	    	     grid_g, & !INTENT(IN) 
	    	     naddsc, & !INTENT(IN) 
	             hw4       !INTENT(IN)	   
  
  USE mem_basic, ONLY: basic_g  !INTENT(IN) 
  
  USE micphys    ,  ONLY: level !INTENT(IN) 
  
  USE rconstants ,  ONLY: cp,p00,cv,rgas,cpi   !INTENT(IN)

  !lfr use mem_aer1, only : aerosol,num_scalar_aer_1st   !INTENT(IN)
  use mem_chem1, only : nspecies_transported        !INTENT(IN)

  use mem_scratch, only  : scratch  ! only scr1, inout
  
  use var_tables, only : scalar_tab & ! (var_p = IN, var_t = INOUT) 
                        ,num_scalar   ! (IN) 
			
  


  !LFR -> To test
  USE stopwatch


  IMPLICIT NONE
  INTEGER,PARAMETER :: advmnt=1

  INTEGER , PARAMETER :: ON=1,OFF=0

  INTEGER , PARAMETER :: use_true_density  = 0 ! 0= OFF, 1=ON
  
  !- for theoretical experiments
  INTEGER , PARAMETER :: theor_wind = 0        ! 0= OFF, 1=ON

  INTEGER           :: mnt_adv_initialized=0
  
  real, parameter :: c1 = cv/rgas, c2 = p00/rgas !c2 = p00*(cpi**c1)/rgas

  TYPE advmnt_vars 
     REAL,POINTER,DIMENSION  (:,:,:)  :: u3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: v3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: w3d 

     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_in 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_out 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_x 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_y 

     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3du 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dv 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dw 

     REAL,POINTER,DIMENSION  (:,:,:)  :: den0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den1_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den2_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den3_3d 

     REAL,POINTER,DIMENSION  (:,:)    :: dxtW 
     REAL,POINTER,DIMENSION  (:,:)    :: dytW 
     REAL,POINTER,DIMENSION  (:)      :: dztW 

  END TYPE advmnt_vars

  TYPE(advmnt_vars), ALLOCATABLE,DIMENSION(:)   :: advmnt_g
  
  INTEGER,PARAMETER :: inorth=1,isouth=2,ieast=3,iwest=4,fromP=1,toP=2
  TYPE messages
     INTEGER :: nPoints
     INTEGER :: iniP
     INTEGER :: endP
  END TYPE messages
  
  TYPE(messages),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: comm
  

  PUBLIC :: advmnt_driver  ! Subroutine
  PUBLIC :: InitCommMatrix ! Subroutine
  PUBLIC :: procMap

  INTEGER, PARAMETER :: nvariables=2

  INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: procMap
  INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:) :: itag1,itag2
  INTEGER, ALLOCATABLE, DIMENSION(:) :: handle

  INTEGER,ALLOCATABLE,DIMENSION(:,:,:) :: west,east,north,south
  INTEGER,ALLOCATABLE,DIMENSION(:,:,:,:,:) :: lineTag
  INTEGER,ALLOCATABLE,DIMENSION(:,:,:,:,:) :: columnTag

  !nvariables=
  INTEGER, PARAMETER :: v_flux=1
  INTEGER, PARAMETER :: v_imxmn=2
  
  
  INTERFACE get
     MODULE procedure &
          get_real, &
          get_int, &
          get_logical
  END INTERFACE

  INTERFACE send
     MODULE procedure &
          send_real, &
          send_int, &
          send_logical
  END INTERFACE
  
  INTEGER :: numberOfProcessors

  type (watchtype) :: w(4)   ! This declares w to be a watch
  type (watchtype) :: w1(6)   ! This declares w to be a watch

CONTAINS
  !----------------------------------------------------
             !call ADVECTc('V' ,mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
  SUBROUTINE advmnt_driver(varn,m1 ,m2 ,m3 ,ia,iz,ja,jz,izu,jzv,mynum)
  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz,ipaths,iget_paths
    IMPLICIT NONE
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    INTEGER , INTENT(IN) :: izu
    INTEGER , INTENT(IN) :: jzv
    INTEGER , INTENT(IN) :: mynum
    character(len=*),intent(IN) :: varn
    
    !--- local vars
    integer n,ng,mxyzp,i,j,procfile
    real, pointer :: scalarp, scalart
    
    INTEGER :: iteste !lfr
    !PRINT *,'No driver da advect ',mynum;CALL flush(6)

    !call create_watch(w,name=(/"Waiting MPI "/))    ! Watches must be created before they are used
    !-scratch arrays initialization  and air density for using
    !- basic state densities
    
    !PRINT *, 'In advection........';CALL flush(6)
    
    iteste=0 !lfr
    
    IF(mnt_adv_initialized == OFF) THEN
       CALL initialize_advmnt(ngrids,nodemzp(mynum,:), &
                              nodemxp(mynum,:),nodemyp(mynum,:))

       CALL prepareCommunication(ia,iz,ja,jz)

       do ng=1,ngrids
         CALL initialize_grid_spacings(ng,nodemzp(mynum,ng), &
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      ,grid_g(ng)%dxt	    &
                      ,grid_g(ng)%dyt	    &
 		      ,grid_g(ng)%fmapt     &
		      ,grid_g(ng)%rtgt      &
!
    		      ,advmnt_g(ng)%dxtW &
    		      ,advmnt_g(ng)%dytW &
		      ,advmnt_g(ng)%dztW )
       enddo
       
       if(use_true_density == OFF) then 
          do ng=1,ngrids
	    CALL initialize_densities(nodemzp(mynum,ng),&
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      , basic_g(ng)%dn0     &
		      , basic_g(ng)%dn0u    &
                      , basic_g(ng)%dn0v    &
                      ,advmnt_g(ng)%dd0_3d  &
    		      ,advmnt_g(ng)%dd0_3du &
    		      ,advmnt_g(ng)%dd0_3dv &
		      ,advmnt_g(ng)%dd0_3dw )
	   enddo
       endif
       
       mnt_adv_initialized= ON    
    ENDIF
    mxyzp=m1*m2*m3
    
    !- Advect  U, V, and W
    if (varn .eq. 'V' .or. varn .eq. 'ALL') then
      stop 'not using mnt to advect u,v,w'
    endif

    !- Advect  scalars
    
    !- get actual air densities, if using them instead of basic state fields
    if(use_true_density == ON) then
     call get_true_densities(m1,m2,m3,level &
                      , basic_g(ngrid)%rtp     &
		      , basic_g(ngrid)%rv      &
                      , basic_g(ngrid)%pp      &
		      , basic_g(ngrid)%pi0     &
                      , basic_g(ngrid)%theta   &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    
    endif

 
    !- prepare wind velocities including map factors
    call prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,basic_g(ngrid)%uc  &
		      ,basic_g(ngrid)%up  &
                      ,basic_g(ngrid)%vc  &
                      ,basic_g(ngrid)%vp  &
		      ,basic_g(ngrid)%wc  &
                      ,basic_g(ngrid)%wp  &
!        
 		      ,grid_g(ngrid)%fmapui &
		      ,grid_g(ngrid)%fmapvi &
		      ,grid_g(ngrid)%rtgt   &
                      ,grid_g(ngrid)%rtgu   &
		      ,grid_g(ngrid)%rtgv   &
		      ,grid_g(ngrid)%f13t   &
                      ,grid_g(ngrid)%f23t   &
!
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  )
                                             
 
    if(theor_wind == on) then
 
        call prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    endif
    !- prepare Walcek's air densities
    
    call get_Walceks_densities(dtlt,m1,m2,m3 &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw &
                      ,advmnt_g(ngrid)%den0_3d &
    		      ,advmnt_g(ngrid)%den1_3d &
    		      ,advmnt_g(ngrid)%den2_3d &
		      ,advmnt_g(ngrid)%den3_3d &
    		      ,advmnt_g(ngrid)%dxtW &
    		      ,advmnt_g(ngrid)%dytW &
		      ,advmnt_g(ngrid)%dztW &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    )

         !call advect_mnt_stalone()
	 
	 !STOP 567

     !- ready to do advection, loop over all scalars
     do n=1,num_scalar(ngrid)     

!srf - somente para gases e aerossois
!     do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
!srf - somente para gases e aerossois
     !do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
        !print*,'nscalar1=',n
	!PRINT *,scalar_tab(n,ngrid)%name
	if (scalar_tab(n,ngrid)%name /= 'COP') cycle
        !print*,'nscalar2=',n
!srf - somente para gases e aerossois
!srf - somente para gases e aerossois
        !PRINT *,n; CALL flush(6)


 !srf-falta incluir controle para shaved_eta e aerossois com sedimentacao


      scalarp => scalar_tab(n,ngrid)%var_p
      scalart => scalar_tab(n,ngrid)%var_t

!for intel - descomente
      call atob(mxyzp,scalarp,advmnt_g(ngrid)%vc3d_in) 
!LFR>       print *,'LFR  (.)(.) atob - passou ', mynum,n,ngrid,num_scalar(ngrid) ; call flush(6)

      call advect_mnt(m1,m2,m3,ia,iz,ja,jz,dtlt &
!for pgi - descomente
!                     ,scalarp                 & ! initial mass mixing ratio (mmxr)
!for intel - descomente
 ,advmnt_g(ngrid)%vc3d_in                     & ! initial mass mixing ratio (mmxr)
 ,advmnt_g(ngrid)%vc3d_out                    & ! final mmxr due advection
 ,advmnt_g(ngrid)%vc3d_x                      &
 ,advmnt_g(ngrid)%vc3d_y		      &
!-- 
                     ,advmnt_g(ngrid)%u3d     & ! 3D U wind
    		     ,advmnt_g(ngrid)%v3d     & ! 3D V wind
    		     ,advmnt_g(ngrid)%w3d     & ! 3D SIGMAZ wind
                     ,advmnt_g(ngrid)%dd0_3d  & ! density 0
                     ,advmnt_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
    		     ,advmnt_g(ngrid)%den1_3d & ! density D1
    		     ,advmnt_g(ngrid)%den2_3d & ! density D2
		     ,advmnt_g(ngrid)%den3_3d & ! density D3
    		     ,advmnt_g(ngrid)%dxtW    &
    		     ,advmnt_g(ngrid)%dytW    &
		     ,advmnt_g(ngrid)%dztW    &
                     ,grid_g(ngrid)%dxt       & ! true dx
                     ,grid_g(ngrid)%dyt       & ! true dy
                     ,dzt(1:m1) 	      & ! true dz
		     ,grid_g(ngrid)%rtgt      &
		     ,mynum,n                    )
     
      call advtndc(m1,m2,m3,ia,iz,ja,jz        &
!                  ,scalarp  ,scratch%scr1 (1)  &
                  ,scalarp  ,advmnt_g(ngrid)%vc3d_out  &
                  ,scalart  ,dtlt,mynum        )

        iteste=iteste+1
     enddo
     
 !call print_watch(w(1))     ! This prints the measured time
      call print_watch(w(1))
     DO i=1,5
        call print_watch(w1(i))
      END DO

      call reset_watch(w) 
      call reset_watch(w1) 

    !PRINT *,'Saindo do driver da advect ',mynum;CALL flush(6)

 
  END SUBROUTINE advmnt_driver
 !----------------------------------------------------
  
  SUBROUTINE initialize_advmnt(ngrids, mmzp,mmxp,mmyp)
   implicit none 
   INTEGER , INTENT(IN) :: ngrids
   INTEGER , INTENT(IN) :: mmxp(ngrids)
   INTEGER , INTENT(IN) :: mmyp(ngrids)
   INTEGER , INTENT(IN) :: mmzp(ngrids)
   
   INTEGER :: ng
   
   !maxz=maxval(mmzp(1:ngrids))
   !maxx=maxval(mmxp(1:ngrids))
   !maxy=maxval(mmyp(1:ngrids))
   
   IF(ALLOCATED(advmnt_g)) THEN 
      PRINT *,'Error in initialize_advmnt, sub: radvc_mnt: advmnt_g already allocated!'
      RETURN
   END IF
   ALLOCATE (advmnt_g(ngrids))

   do ng=1,ngrids
      ALLOCATE(advmnt_g(ng)%u3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%u3d=0.
      ALLOCATE(advmnt_g(ng)%v3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%v3d=0.
      ALLOCATE(advmnt_g(ng)%w3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%w3d=0.
      
      ALLOCATE(advmnt_g(ng)%dd0_3d (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3d =0.
      ALLOCATE(advmnt_g(ng)%dd0_3du(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3du=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dv(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dv=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dw(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmnt_g(ng)%den0_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den0_3d=0.
      ALLOCATE(advmnt_g(ng)%den1_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den1_3d=0.
      ALLOCATE(advmnt_g(ng)%den2_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den2_3d=0.
      ALLOCATE(advmnt_g(ng)%den3_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den3_3d=0.

      ALLOCATE(advmnt_g(ng)%dxtW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dxtW=0.
      ALLOCATE(advmnt_g(ng)%dytW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dytW=0.
      ALLOCATE(advmnt_g(ng)%dztW(mmzp(ng)))                  ; advmnt_g(ng)%dztW=0.


      ALLOCATE(advmnt_g(ng)%vc3d_in  (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_in =0.
      ALLOCATE(advmnt_g(ng)%vc3d_out (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_out=0.
      ALLOCATE(advmnt_g(ng)%vc3d_x   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_x  =0.
      ALLOCATE(advmnt_g(ng)%vc3d_y   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_y  =0.


   enddo

   



  END SUBROUTINE initialize_advmnt
 !----------------------------------------------------
 
  SUBROUTINE initialize_densities(m1,m2,m3,dn0,dn0u,dn0v &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::dn0,dn0u,dn0v
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::dd0_3d,dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   
   dd0_3d (:,:,:)=  dn0 (:,:,:) 
   dd0_3du(:,:,:)=  dn0u(:,:,:)
   dd0_3dv(:,:,:)=  dn0v(:,:,:)
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dn0(k,i,j) +dn0(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo
   

  END SUBROUTINE initialize_densities
 !----------------------------------------------------
 
  SUBROUTINE initialize_grid_spacings(ng,m1,m2,m3,dxt,dyt,fmapt,rtgt &
			    ,dxtW,dytW,dztW )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: ng,m1,m2,m3
   REAL,DIMENSION(m2,m3),intent(IN) :: dxt,dyt,fmapt,rtgt
   
   REAL,DIMENSION(m2,m3),intent(OUT):: dxtW,dytW
   REAL,DIMENSION(m1),intent(OUT):: dztW
   ! local var
   integer i,j,k
   real rtgti
   
   do j = 1,m3
      do i = 1,m2
   	rtgti = 1. / rtgt(i,j)
        
	!- at init/rams_grid.f90:
        !     dxt(i,j)=fmapt(i,j)/(xmn(i,ngrid)-xmn(i-1,ngrid))  
        !     dyt(i,j)=fmapt(i,j)/(ymn(j,ngrid)-ymn(j-1,ngrid))

   	 dxtW(i,j) = 1./(dxt(i,j) * fmapt(i,j) * rtgti)
   	 dytW(i,j) = 1./(dyt(i,j) * fmapt(i,j) * rtgti)
   enddo;enddo
   do k = 1,m1
    !- at init/gridset.f90:	
    !  dztn(k,ifm) = 1. / (zmn(k,ifm) - zmn(k-1,ifm))
    ! Por que o Jacobiano nao depende de Z, o dztw depende somente 
    ! de z.
    !dztW(k,i,j) = 1./ ( dzt(k) * rtgti * fmapt(i,j)**2 )
     dztW(k)	 = 1./ ( dztn(k,ng) ) !
 
   enddo
  END SUBROUTINE initialize_grid_spacings
 !----------------------------------------------------
 
  SUBROUTINE get_true_densities(m1,m2,m3,level,rtp,rv,pp,pi0,theta &
			       ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,level
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::rtp,rv,pp,pi0,theta
   REAL,DIMENSION(:,:,:),intent(OUT):: dd0_3d
   REAL,DIMENSION(m1,m2,m3),intent(OUT)  :: dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   real c3
   
   c3 = c2 * (cpi**c1)

   !- true air density at points "T"
   
   if( level == 0 ) then
     dd0_3d(:,:,:) = (c3/theta(:,:,:))*(pi0(:,:,:)+pp(:,:,:))**c1
   else
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1
            dd0_3d(k,i,j) = (c3/theta(k,i,j))* (1. + rtp(k,i,j))/ &
	            (1. + 1.61*rv(k,i,j))*(pi0(k,i,j)+pp(k,i,j))**c1
         end do
      end do
   end do
   endif 

   !- true air density at points "U", "V" and "W":
   
   call fill_dn0uv(m1,m2,m3,dd0_3d,dd0_3du,dd0_3dv)

   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dd0_3d(k,i,j) + dd0_3d(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo

 END SUBROUTINE get_true_densities

 !----------------------------------------------------

 SUBROUTINE prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,uc,up,vc,vp,wc,wp &
  		            ,fmapui &
		            ,fmapvi &
                            ,rtgt   &
                            ,rtgu   &
		            ,rtgv   &
		            ,f13t   &
                            ,f23t   &
			    ,u3d,v3d,w3d)
   
   implicit none 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL,DIMENSION(m1,m2,m3),intent(INOUT) :: uc,up,vc,vp,wc,wp
!   REAL,DIMENSION(m1,m2,m3),intent(IN) :: uc,up,vc,vp,wc,wp ! TMP <<<<<<<<<<<<<<<<<<<<<<<
   REAL,DIMENSION(m2,m3)   ,intent(IN) :: rtgt,rtgu,rtgv,fmapui,fmapvi,f13t,f23t
   
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d
  
   !- local var
   real   dtlto2
   integer jm,jp,im,ip 
   integer i,j,k
   real :: c1,c2,rtgti

   dtlto2 = .5! * dtlt

  ! u3d, u3d, and w3d are input as the velocity components (averaged
  ! between past and current time levels) times dtlt.
   do j=1,m3
     do i = 1,m2
      do k = 1,m1
  
          w3d(k,i,j) = ( wc(k,i,j) + wp(k,i,j) )*dtlto2
          u3d(k,i,j) = ( uc(k,i,j) + up(k,i,j) )*dtlto2
          v3d(k,i,j) = ( vc(k,i,j) + vp(k,i,j) )*dtlto2
  
   enddo;enddo;enddo


  !return ! for pure cartesian coordinates

  ! here w3d is the cartesian vertical velocity 
  
  ! Add contribution to w3d from horiz winds crossing sloping sigma surfaces,
  ! and include 1/rtgt factor in w3d
  do j = 1,m3		 
     jm = max(1,j-1)
     jp = min(m3,j+1)
     do i = 1,m2
  	im = max(1,i-1)
  	ip = min(m2,i+1)
        rtgti = 1. / rtgt(i,j)
	
	do k = 1,m1-1
	    w3d(k,i,j) = ( (u3d(k,i,j) + u3d(k+1,i,j) + u3d(k,im,j) + u3d(k+1,im,j) ) * f13t(i,j)  &
	               +   (v3d(k,i,j) + v3d(k+1,i,j) + v3d(k,i,jm) + v3d(k+1,i,jm) ) * f23t(i,j)  &
		         ) * hw4(k)  &
	               + w3d(k,i,j) * rtgti

  	enddo
     enddo
  enddo
  ! here w3d is the sigma_z velocity 

  !- including map factors on U,V:
  do j = 1,m3
     do i = 1,m2
  	c1 = fmapui(i,j) * rtgu(i,j)
  	c2 = fmapvi(i,j) * rtgv(i,j)
  	do k = 1,m1-1
  	   u3d(k,i,j) = u3d(k,i,j) * c1 
  	   v3d(k,i,j) = v3d(k,i,j) * c2 
	   
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	   !u3d(k,i,j) = 0. !<<<<<<<<<<<<<<<<<<<<<<
	   
	   
	   !v3d(k,i,j) = -1. !<<<<<<<<<<<<<<<<<<<<<<

           !w3d(k,i,j) =	 0.
	   !wc(k,i,j)=w3d(k,i,j)
	   !uc(k,i,j)=u3d(k,i,j)
	   !vc(k,i,j)=v3d(k,i,j)
	   !wp(k,i,j)=w3d(k,i,j)
	   !up(k,i,j)=u3d(k,i,j)
	   !vp(k,i,j)=v3d(k,i,j)
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  	enddo
     enddo
  enddo

  END SUBROUTINE prepare_winds
 !----------------------------------------------------

  SUBROUTINE get_Walceks_densities(dt,m1,m2,m3,u3d,v3d,w3d &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw &
			    ,den0_3d,den1_3d,den2_3d,den3_3d &
			    ,dxtW,dytW,dztW,dxt,dyt)

   IMPLICIT NONE 
   !-in
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL    , INTENT(IN) :: dt
   REAL, DIMENSION(m1)      , INTENT(IN) :: dztW
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL, DIMENSION(m1,m2,m3), INTENT(IN) :: u3d,v3d,w3d     &
			                   ,dd0_3du &
					   ,dd0_3dv,dd0_3dw
   REAL, DIMENSION(m1,m2,m3), INTENT(INOUT) :: dd0_3d 
   !-out
   REAL, DIMENSION(m1,m2,m3), INTENT(OUT):: den0_3d,den1_3d &
                                           ,den2_3d,den3_3d
   
   ! local var
   integer i,j,k
   
    DO  j=m3,2,-1
      DO  i=2,m2
        DO k = 2,m1
	

!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!dd0_3d(k,i,j)=1.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


            den0_3d(k,i,j)=dd0_3d(k,i,j)



	    
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxt(i,j)*(d0         *u3d(k,i,j)- d0           *u3d(k,i-1,j ))
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt*dxt(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))
	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxtW(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))

!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dx*(d0*v3d(k,i,j)- d0*v3d(k,i,j-1 ))
!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt*dyt(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
            den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dytW(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
!
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dx         *(d0*w3d(k,i,j)-d0*w3d(  k-1,i,j))
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt*dzt(k)     *(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))
            den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dztW(k)*(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))

       END DO !i
      END DO !j
    END DO !k  

 

 END SUBROUTINE get_Walceks_densities
 !----------------------------------------------------
 SUBROUTINE advect_mnt(m1,m2,m3,ia,iz,ja,jz,dt,vc3d_in,vc3d_out,vc3d_x,vc3d_y&  
                      ,u3d,v3d,w3d,dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d &
		      ,dxtW,dytW,dztW,dxt,dyt,dzt,rtgt,mynum,n)
   
   use node_mod, only : nodei0,nodej0
  
  IMPLICIT NONE
  INTEGER , INTENT(IN) :: m1
  INTEGER , INTENT(IN) :: m2
  INTEGER , INTENT(IN) :: m3
  INTEGER , INTENT(IN) :: ia
  INTEGER , INTENT(IN) :: iz
  INTEGER , INTENT(IN) :: ja
  INTEGER , INTENT(IN) :: jz,n
  INTEGER , INTENT(IN) :: mynum
  
  REAL    , INTENT(IN) :: dt
  
  REAL,DIMENSION(m1)      , INTENT(IN)  :: dzt
  REAL,DIMENSION(m2,m3)   , INTENT(IN)  :: dxt,dyt,rtgt
  
  REAL,DIMENSION(m1),       INTENT(IN) :: dztW 
  REAL,DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

  REAL,DIMENSION(m1,m2,m3), INTENT(INOUT) :: vc3d_in! initial mixing ratio 
  REAL,DIMENSION(m1,m2,m3), INTENT(INOUT) :: u3d,v3d,w3d
  REAL,DIMENSION(m1,m2,m3), intent(INOUT) :: dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d

  REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_x,vc3d_y
  REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_out ! updated mixing ratio due advection

  !- local var
  REAL,DIMENSION(m1)               :: dxx  
  REAL,DIMENSION(m2,m3)            :: dxy  
  real masscon,initialmass,vol
  integer :: i,j,k,nrec,it
  INTEGER,SAVE :: step=0
  
  CHARACTER(LEN=2) :: nf

  STEP=STEP+1

!  !For test standalone write
!  WRITE(nf,FMT='(I2.2)') step
!  OPEN(27,file='walcekDataIn.'//nf,form='unformatted',STATUS='replace')
!  WRITE(27) m1,m2,m3,ia,iz,ja,jz
!  WRITE(27) dt,nodei0,nodej0
!  WRITE(27) dzt,dxt,dyt,rtgt,dztW,dxtW,dytW
!  WRITE(27) vc3d_in,u3d,v3d,w3d,dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d
!  PRINT *,'Done!'
!  CLOSE(27)
!  !End Test write state
  

!---------------------------------------------------------------------  
!- for theoretical experiments
!
!  OPEN(19,file='walcek_mixingratio_3d_rot_wind.gra',         &
!       form='unformatted',access='direct',status='unknown',  &
!       recl=4*m2*m3)
!  nrec =0
!  do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) u3d(k,:,:)
!  enddo
!  do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) v3d(k,:,:)
!  enddo
!  CLOSE(19)
!
!  OPEN(19,file='walcek_mixingratio_3d_rot.gra',         &
!       form='unformatted',access='direct',status='unknown',  &
!       recl=4*m2*m3)
!  nrec =0
!  
!  DO  it=1,3500
!dd0_3d
!   do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) vc3din(k,:,:)
!   enddo
! !-------------------------------
!  !- check mass conservation (initial mass)
!  initialmass=0.
!  do j=2,m3-1; do i=2,m2-1
!    do k=2,m1-1
!    !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
!    vol= 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j) ! true volume
!    initialmass=initialmass+vc3din(k,i,j)*dd0_3d(k,i,j)*vol
!
!  enddo;enddo;enddo
!---------------------------------------------------------------------
!---------------------------------------------------------------------


  !vc3d_x=vc3d_in
  !vc3d_y=vc3d_in
  !vc3d_out=vc3d_in
  
  
  ! Do I-advection first over all J-k lines
 
  !WRITE (80+mynum,FMT='(A)') 'LFR (.) -> X';CALL flush(80+mynum)

  !call Advec3d_X(m1,m2,m3,ia,iz,ja,jz,vc3din,u3d,den0_3d, den1_3d,dt,dxtW,dd0_3d,vc3din)
  !print *, 'Advect em x ',mynum;CALL flush(6)
  !go to 500
  call start_watch(w(4))
  call start_watch(w(1))
  call Advec3d_X(m1,m2,m3,ia,iz,ja,jz,vc3d_in,u3d,den0_3d, &
                 den1_3d,dt,dxtW,dd0_3d,vc3d_x,mynum,n )
  call stop_watch(w(1))
  !srf - only X-adv
  !vc3d_out=vc3d_x
  !return
  
  !500 continue
  
  !- do J-advection next over all I-k lines
  !call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3din,v3d,den1_3d, den2_3d,dt,dytW,dd0_3d,vc3din)
  !print *, 'Advect em y ',mynum;CALL flush(6)
 
  !vc3d_out=vc3d_in
  !vc3d_x=vc3d_in ! <<<<<< TMP
  !call start_watch(w(2))
  call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3d_x ,v3d,den1_3d, &
                 den2_3d,dt,dytW,dd0_3d,vc3d_y )
  !call stop_watch(w(2))
  !srf - only Y-adv
  !vc3d_out=vc3d_y
  !return
  
  
  !- do k-advection next over all I-J lines
  ! vc3dout_Z = razao de mistura final (apos a adveccao nas 3 dimensoes)
  !call Advec3d_Z(m1,m2,m3,ia,iz,ja,jz,vc3din,w3d,den2_3d, den3_3d,dt,dxx,dd0_3d,vc3dout_Z)
  !print *, 'Advect em z ',mynum;CALL flush(6)
  !call start_watch(w(3))
  call Advec3d_Z(m1,m2,m3,ia,iz,ja,jz,vc3d_y,w3d,den2_3d, &
                 den3_3d,dt,dztW,dd0_3d,vc3d_out)
  
  
!   !For test standalone
!  OPEN(19,file='walcekDataOut.'//nf//'.gra',         &
!       form='unformatted',access='direct',status='unknown',  &
!       recl=4*m2*m3)
!  nrec =0
!  do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) vc3d_out(k,:,:)
!  enddo
!  CLOSE(19)
!  IF (step==20) STOP 455 
! !End test

  
  
  !call stop_watch(w(3))
  call stop_watch(w(4))
  !WRITE (80+mynum,FMT='(A)') 'LFR (.) -> Fim';CALL flush(80+mynum)
 
  !-------------------------------
  !- check mass conservation (after advection)
  !masscon=0.
  !do j=2,m3-1; do i=2,m2-1
  !  do k=2,m1-1
  !  !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
  !  vol= 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j) ! true volume
  !  masscon=masscon+vc3dout_Z(k,i,j)*dd0_3d(k,i,j)*vol
  !
  !enddo;enddo;enddo
  !print*,'masscon=',initialmass,masscon,100.*(masscon-initialmass)/initialmass
  !print*,'MAXVALS=', maxval(vc3din)*28.96/28, maxval(vc3dout_Z)*28.96/28
  !call flush(6)
  !-------------------------------

!---------------------------------------------------------------------
!- for theoretical experiments
!---------------------------------------------------------------------
!  vc3din=vc3dout_Z
!
! ENDDO
! close (19)
! stop 333
!---------------------------------------------------------------------
!---------------------------------------------------------------------
   
 END SUBROUTINE advect_mnt
 !----------------------------------------------------
 
SUBROUTINE Advec3d_X(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn,mynum,nvar)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !

  use node_mod, only : nodei0,nodej0
 
  IMPLICIT none  
 
  INCLUDE "i8.h"
 
  INTEGER,INTENT(IN)                        :: m1,m2,m3, ia,iz,ja,jz,mynum,nvar
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: den1
  REAL   ,INTENT(in)                        :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)       :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)   :: qn

  REAL,DIMENSION(m1,m2,m3)    :: flux
  REAL,DIMENSION(m1,m2,m3)    :: vcmax
  REAL,DIMENSION(m1,m2,m3)    :: vcmin
  
  LOGICAL,DIMENSION(m1,m2,m3) :: imxmn
  
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: idime
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n,cs

  integer, parameter :: debug =0, dump=0
  integer, parameter :: on = 1, up=1,down=0
  
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz
  
  integer :: nf,itag
  

  REAL, DIMENSION(2*m1,ja:jz) :: buffer
    call start_watch(w1(5))
  nf=80+mynum
  idime = m2 ! x dir

  qn=q0

  imxmn=.false.
  flux=0.0
  ! imxmn(idime+1)=.false.
  
  !WRITE (nf,FMT='(A)') 'X Begin'; CALL flush(nf)
  
  ii=ia+nodei0(mynum,ngrid)
  ii0=ia+nodei0(mynum,ngrid)-1
  ie=iz+nodei0(mynum,ngrid)
  ie0=iz+nodei0(mynum,ngrid)+1
  !print *,nvar
    call stop_watch(w1(5))
    
  call start_watch(w1(1)) !waiting time up
  DO j=ja,jz

    DO k=2,m1-1 
       ! Update mixing ratios and limit Fluxes going UP where u>0
       !  First assume upstream flux at edge of domain
       IF(u(k,1,j)>=zr0) flux(k,1,j)= q0(k,1,j)*u(k,1,j)*dt*dd0(k,1,j)
    END DO

    ji=j+nodej0(mynum,ngrid)
    IF(north(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,1)
        IF (dump==on) WRITE (nf,FMT='("Get - Np: ",I2.2," J: ",I3.3," North: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,north(ngrid,mynum,ji),itag;CALL flush(nf)
        CALL get(buffer(:,j),2*m1,north(ngrid,mynum,ji),itag)
     END IF
   END DO 
   call stop_watch(w1(1))  

   DO j=ja,jz
    call start_watch(w1(1))
    ji=j+nodej0(mynum,ngrid)
    IF(north(ngrid,mynum,ji)/=0) THEN
       itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,1)
       IF (dump==on)  WRITE (nf,FMT='("Wait - Np: ",I2.2," J: ",I3.3," Tag: ",I5.5)') &
	                     mynum,j,itag;CALL flush(nf)
       CALL parf_wait_nostatus(handle(itag))
       !IF (dump==on) WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
       CALL unPackBuffer(m1,buffer(:,j),flux(:,1,j),imxmn(:,1,j))
       IF (dump==on) CALL dumpData(m1,nf,flux(:,1,j),imxmn(:,1,j),buffer(:,j))
    END IF
    call stop_watch(w1(1))
   ! Identify local max and min, specify mixing ratio limits at new time
    call start_watch(w1(3))
  DO  i=2,idime-1 ! ia,iz-1 or 1,iz-1

     DO k=2,m1-1 
				  
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i-1,j),q0(k,i+1,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i-1,j),q0(k,i+1,j))       !       extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i+1,j)
        if(u(k,i-1,j)>=zr0) ck2= q0(k,i-1,j)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !     mixing ratio at t+dt. If these limits are ever violated,
        !     non-monotonic (oscillatory) behavior in solution results
        !

        IF(u(k,i,j)<zr0) CYCLE
     
        IF(u(k,i-1,j)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i+1,j)-q0(k,i-1,j))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i+1,j))  ), MAX(q0(k,i,j),q0(k,i+1,j)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(vcmin(k,i,j),MIN(vcmax(k,i,j),          &   !eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,I-1,j)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i-1,j)
        END IF                                                  !Eq-9a
     END DO
     
  END DO
  call stop_watch(w1(3))
END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  call start_watch(w1(4))
  DO j=ja,jz

     ji=j+nodej0(mynum,ngrid)
     IF(south(ngrid,mynum,ji)/=0) THEN
         itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,2)
         CALL packBuffer(m1,flux(:,idime-1,j),imxmn(:,idime-1,j),buffer(:,j))
         IF (dump==on) WRITE (nf,FMT='("Put - Np: ",I2.2," J: ",I3.3," South: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,south(ngrid,mynum,ji),itag;CALL flush(nf)	 
	 IF (dump==on) CALL dumpData(m1,nf,flux(:,idime-1,j),imxmn(:,idime-1,j),buffer(:,j))
        CALL send(buffer(:,j),2*m1,south(ngrid,mynum,ji),itag)
     END IF  
  END DO
  call stop_watch(w1(4))
   !IF (nvar==3) THEN
   !    WRITE (nf,FMT='(A,2(I3.3,1X))')  'flux(1:2,:,7) &  imxmn(1:2,:,7)  --- ', &
    !                                    1+nodei0(mynum,ngrid),m2+nodei0(mynum,ngrid)
    !   DO k=1,2
!	  WRITE (nf,FMT='(I3.3,1X,$)') k!
!	  DO i=1,m2
!	     WRITE (nf,FMT='(E15.6,1X,L1,1X,$)') flux(k,i,7),imxmn(k,i,7)
 !         END DO
!	  WRITE (nf,FMT='A') ''
!       END DO
!       CALL flush(nf)
!    END IF

 DO j=ja,jz

     call start_watch(w1(2))
     ji=j+nodej0(mynum,ngrid)
     IF(south(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,1)
        IF (dump==on) WRITE (nf,FMT='("Get - Np: ",I2.2," J: ",I3.3," South: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,south(ngrid,mynum,ji),itag;CALL flush(nf)
        CALL get(buffer(:,j),2*m1,south(ngrid,mynum,ji),itag)
     ELSE
        DO k=2,m1-1       
           IF(u(k,idime-1,j)<zr0) flux(k,idime-1,j)= &
    	       q0(k,idime,j)*u(k,idime-1,j)*dt*dd0(k,idime-1,j)
        END DO
     END IF
  END DO 
  call stop_watch(w1(2))
  DO j=ja,jz
    ji=j+nodej0(mynum,ngrid)
    call start_watch(w1(2))
    IF(south(ngrid,mynum,ji)/=0) THEN
	itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,1)
          IF (dump==on) WRITE (nf,FMT='("Wait - Np: ",I2.2," J: ",I3.3," Tag: ",I5.5)') &
	                     mynum,j,itag;CALL flush(nf)
	  CALL parf_wait_nostatus(handle(itag))
          !WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
          CALL unPackBuffer(m1,buffer(:,j),flux(:,idime-1,j),imxmn(:,idime,j))
	  IF (dump==on) CALL dumpData(m1,nf,flux(:,idime-1,j),imxmn(:,idime,j),buffer(:,j))
     END IF
     call stop_watch(w1(2))
     call start_watch(w1(3))
     DO i=idime-1,2,-1 !iz,ia,-1
       ! print*,'part 2'
       DO k=2,m1-1 
        
        IF(u(k,i-1,j)>=zr0) THEN           ! Inflow-only cell
        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))
	ELSE
              x1=  dt*ABS(u(k,i-1,j))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i-1,j)-q0(k,i+1,j))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i-1,j)) ), MAX(q0(k,i,j),q0(k,i-1,j)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i-1,j))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,I,j)/DXX(I,J)-X1*CF1*DD0(k,I-1,j))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i-1,j)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
           
	     !write(mynum,*) i+nodei0(mynum,ngrid),qn(k,i,j),q0(k,i,j),flux(k,i,j),flux(k,i-1,j); call flush(mynum)
             !if(i+nodei0(mynum,ngrid) == 16) then
   	     !  write(mynum,*)  vcmin(k,i,j), vcmax(k,i,j),q0(k,i,j),cf
   	     ! 	write(mynum,*)  x1,cf1,dd0(k,i-1,j),den0(k,i,j),den1(k,i,j)
             ! endif     ; call flush(mynum)
        END IF                                                  
     END DO

  END DO
      call stop_watch(w1(3))
  !END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop y-z
1000 continue
 call start_watch(w1(4))

  DO j=ja,jz

    je=j+nodej0(mynum,ngrid)
    ji=j+nodej0(mynum,ngrid)

     !sending last column  (i=idime-1)  to east node
     IF(north(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,2)
        CALL packBuffer(m1,flux(:,1,j),imxmn(:,2,j),buffer(:,j))
	IF (dump==on) WRITE (nf,FMT='("Put - Np: ",I2.2," J: ",I3.3," North: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,North(ngrid,mynum,ji),itag;CALL flush(nf)	
        CALL send(buffer(:,j),2*m1,north(ngrid,mynum,ji),itag)
	IF (dump==on) CALL dumpData(m1,nf,flux(:,1,j),imxmn(:,2,j),buffer(:,j))
     END IF  
  END DO
   call stop_watch(w1(4))
  !WRITE (nf,FMT='(A)') 'X ended'; CALL flush(nf)
  !STOP 500
END SUBROUTINE Advec3d_X
!--------------------------------------------------------------------------

SUBROUTINE Advec3d_Y(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  use node_mod, only : nodei0,nodej0

  IMPLICIT none  
   
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                       :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1,m2,m3)                 :: flux
  REAL,DIMENSION(m1,m2,m3)                   :: vcmax
  REAL,DIMENSION(m1,m2,m3)                   :: vcmin
  LOGICAL,DIMENSION(m1,m2,m3)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n,cs

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz,nvar
  
  
  integer :: nf,itag

  REAL, DIMENSION(2*m1,ia:iz) :: buffer
  nf=80+mynum

  idime = m3 ! y dir

  qn= q0
  imxmn=.false.
  flux=0.0
  !WRITE (nf,FMT='(A)') 'Y Begin'; CALL flush(nf)
  
  ji =ja+nodej0(mynum,ngrid)
  ji0=ja+nodej0(mynum,ngrid)-1
  je =jz+nodej0(mynum,ngrid)
  je0=jz+nodej0(mynum,ngrid)+1

  DO i=ia,iz

    ! Update mixing ratios and limit Fluxes going UP where u>0
    !  First assume upstream flux at edge of domain
    DO k=2,m1-1   
       IF(u(k,i,1)>=zr0) flux(k,i,1)= q0(k,i,1)*u(k,i,1)*dt*dd0(k,i,1)
    END DO
      
    ii=i+nodei0(mynum,ngrid)

    IF(west(ngrid,mynum,ii)/=0) THEN
       itag=lineTag(ngrid,mynum,west(ngrid,mynum,ii),ii,1)
       !WRITE (nf,FMT='("Get - Np: ",I2.2," I: ",I3.3," West : ",I2.2," Tag: ",I5.5)') &
	!                     mynum,i,west(ngrid,mynum,ii),itag;CALL flush(nf)
        CALL get(buffer(:,i),2*m1,west(ngrid,mynum,ii),itag)
    END IF

  END DO 

 !- big loop y-z
 DO i=ia,iz
    ii=i+nodei0(mynum,ngrid)

     IF(west(ngrid,mynum,ii)/=0) THEN
       itag=lineTag(ngrid,mynum,west(ngrid,mynum,ii),ii,1)
       !WRITE (nf,FMT='("Wait - Np: ",I2.2," I: ",I3.3," Tag: ",I5.5)') &
	!                     mynum,i,itag;CALL flush(nf)
        CALL parf_wait_nostatus(handle(itag))
        !WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
        CALL unPackBuffer(m1,buffer(:,i),flux(:,i,1),imxmn(:,i,1))	   
     END IF
  
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  j=2,idime-1 ! ja,jz
     DO k=2,m1-1   
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i,j-1),q0(k,i,j+1)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i,j-1),q0(k,i,j+1))	    !	    extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i,j+1)
        if(u(k,i,j-1)>=zr0) ck2= q0(k,i,j-1)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
!     END DO
!  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  
!  DO j=2,idime-1 ! ja,jz
!     DO k=2,m1-1     
        !print*,'part 1'
        IF(u(k,i  ,j)<zr0) CYCLE
     
        IF(u(k,i,j-1)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i,j+1)-q0(k,i,j-1))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j+1))  ), MAX(q0(k,i,j),q0(k,i,j+1)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),          &   !eq-3&8
                 (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,i,J-1)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	     !write(mynum,*) j+nodej0(mynum,ngrid),qn(k,i,j),q0(k,i,j),flux(k,i,j),flux(k,i,j-1); call flush(mynum)
             !if(j+nodej0(mynum,ngrid) == 30) then
   	     !  write(mynum,*)  vcmin(k,i,j), vcmax(k,i,j),q0(k,i,j),cf
   	     !	write(mynum,*)  x1,cf1,dd0(k,i,j-1),den0(k,i,j),den1(k,i,j)
             !endif     ; call flush(mynum)
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i,j-1)
        END IF                                                  !Eq-9a
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  END DO
  DO i=ia,iz
     ii=i+nodei0(mynum,ngrid)
     IF(east(ngrid,mynum,ii)/=0) THEN
        itag=lineTag(ngrid,mynum,east(ngrid,mynum,ii),ii,2)
	CALL packBuffer(m1,flux(:,i,idime-1),imxmn(:,i,idime-1),buffer(:,i))
        !WRITE (nf,FMT='("Put - Np: ",I2.2," i: ",I3.3," east: ",I2.2," Tag: ",I5.5)') &
	!                     mynum,i,east(ngrid,mynum,ii),itag;CALL flush(nf)
	CALL send(buffer(:,i),2*m1,east(ngrid,mynum,ii),itag)
     END IF  
  END DO
  
 DO i=ia,iz

    ii=i+nodei0(mynum,ngrid)

    IF(east(ngrid,mynum,ii)/=0) THEN
	itag=lineTag(ngrid,mynum,east(ngrid,mynum,ii),ii,1)
        !WRITE (nf,FMT='("Get - Np: ",I2.2," I: ",I3.3," east : ",I2.2," Tag: ",I5.5)') &
	!                     mynum,i,east(ngrid,mynum,ii),itag;CALL flush(nf)
        CALL get(buffer(:,i),2*m1,east(ngrid,mynum,ii),itag)
     ELSE
       DO k=2,2!2,m1-1   
          IF(u(k,i,idime-1)<zr0) flux(k,i,idime-1)=q0(k,i,idime)*u(k,i,idime-1)*dt*dd0(k,i,idime-1)
       END DO

     END IF

  END DO 

  DO i=ia,iz
     ii=i+nodei0(mynum,ngrid)
     IF(east(ngrid,mynum,ii)/=0) THEN
	itag=lineTag(ngrid,mynum,east(ngrid,mynum,ii),ii,1)     
        !WRITE (nf,FMT='("Wait - Np: ",I2.2," I: ",I3.3," Tag: ",I5.5)') &
	!                     mynum,i,itag;CALL flush(nf)	  
	CALL parf_wait_nostatus(handle(itag))
	!WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
	CALL unPackBuffer(m1,buffer(:,i),flux(:,i,idime-1), &
	                 imxmn(:,i,idime))
     END IF

  
  DO j=idime-1,2,-1 !jz,ja,-1
     DO k=2,m1-1     
        IF(u(k,i,j-1)>=zr0) THEN           ! Inflow-only cell        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))
        ELSE
              x1=  dt*ABS(u(k,i,j-1))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i,j-1)-q0(k,i,j+1))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j-1)) ), MAX(q0(k,i,j),q0(k,i,j-1)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                     (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i,j-1))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,i,J)/DXX(I,J)-X1*CF1*DD0(k,I,j-1))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i,j-1)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
	     
	     !write(mynum,*) j+nodej0(mynum,ngrid),qn(k,i,j),q0(k,i,j),flux(k,i,j),flux(k,i,j-1); call flush(mynum)
             !if(j+nodej0(mynum,ngrid) == 15) then
   	     !  write(mynum,*)  vcmin(k,i,j), vcmax(k,i,j),q0(k,i,j),cf
   	     !	write(mynum,*)  x1,cf1,dd0(k,i,j-1),den0(k,i,j),den1(k,i,j)
             !endif     ; call flush(mynum)
        END IF                                                  
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop x-z

  DO i=ia,iz
     ii=i+nodei0(mynum,ngrid)
     IF(west(ngrid,mynum,ii)/=0) THEN
	itag=lineTag(ngrid,mynum,west(ngrid,mynum,ii),ii,2)
        CALL packBuffer(m1,flux(:,i,1),imxmn(:,i,2),buffer(:,i))
        !WRITE (nf,FMT='("Put - Np: ",I2.2," i: ",I3.3," West: ",I2.2," Tag: ",I5.5)') &
	!                     mynum,i,west(ngrid,mynum,ii),itag;CALL flush(nf)
	CALL send(buffer(:,i),2*m1,west(ngrid,mynum,ii),itag)
     END IF  
  END DO
  !WRITE (nf,FMT='(A)') 'Y ended'; CALL flush(nf)

END SUBROUTINE Advec3d_Y
!---------------------------------------------------------------------------



SUBROUTINE Advec3d_Z(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m1)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1)                 :: flux
  REAL,DIMENSION(m1)                   :: vcmax
  REAL,DIMENSION(m1)                   :: vcmin
  LOGICAL,DIMENSION(m1)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

 idime = m1 ! z dir

  qn = q0
  !PRINT *, 'LFR -> ia,iz,m2,ja,jz,m3,m1=',ia,iz,m2,ja,jz,m3,m1;CALL flush(6)
 !- big loop y-x
 DO j=ja,jz ; DO i=ia,iz
 
  
  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  k=2,idime-1 !
  
     imxmn(k)=q0(k,i,j)>=max(q0(k-1,i,j),q0(k+1,i,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k-1,i,j),q0(k+1,i,j))	    !	    extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k  ,i,j)< zr0) ck1= q0(k+1,i,j)
     if(u(k-1,i,j)>=zr0) ck2= q0(k-1,i,j)
     vcmax(k)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(k)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1,i,j)>=zr0) flux(1)= q0(1,i,j)*u(1,i,j)*dt*dd0(1,i,j)
  
  DO k=2,idime-1 ! ja,jz
  
  !print*,'part 1'
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k-1,i,j)<zr0) THEN
        flux(k)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream
        x1= dt*u(k,i,j)/dxx(k)               ! Courant number
        x1n= (1.-x1)*(q0(k+1,i,j)-q0(k-1,i,j))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(k-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(k+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k+1,i,j))  ), MAX(q0(k,i,j),q0(k+1,i,j)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(k),   MIN(   vcmax(k),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k-1)/dxx(k))/den1(k,i,j) ))

	! for debug purposes only
        if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k-1)/dxx(k))/DEN1(k,i,j)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(k)= dxx(k)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  
  IF(u(idime-1,i,j)<zr0) flux(idime-1)=q0(idime,i,j)*u(idime-1,i,j)*dt*dd0(idime-1,i,j)
 
  DO k=idime-1,2,-1 !jz,ja,-1
     
       IF(u(k-1,i,j)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k),   MIN(   vcmax(k),&
           (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k) + flux(k-1)/dxx(k))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k-1,i,j))/dxx(k)     ! Courant number
           x1n= (1.-x1)*(q0(k-1,i,j)-q0(k+1,i,j))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(k+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k-1,i,j)) ), MAX(q0(k,i,j),q0(k-1,i,j)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(k),  MIN(   vcmax(k), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k)-x1*cf1*dd0(k-1,i,j))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k)/dxx(k)-X1*CF1*DD0(k-1,I,j))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(k-1)=dxx(k)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop y-x


END SUBROUTINE Advec3d_Z

  SUBROUTINE dumpData(m1,nf,var1,var2,var3)
    INTEGER, INTENT(IN) :: m1,nf
    REAL, INTENT(IN) :: var1(m1)
    LOGICAL,INTENT(IN) :: var2(m1)
    REAL, INTENT(IN) :: var3(2*m1)
   
    INTEGER :: k
    REAL :: cs1,cs2
   
    cs1=0.0;cs2=0.0
   
    DO k=2,m1-1; WRITE (nf,FMT='(E15.5,1X,$)') var1(k); cs1=cs1+var1(k) ; END DO; WRITE (nf,FMT='(A)') ''; CALL flush(nf)
    DO k=2,m1-1
       WRITE (nf,FMT='(L1,1X,$)') var2(k)
       IF(var2(k)) THEN
          cs2=cs2+1
       END IF 
    END DO
    WRITE (nf,FMT='(A)') ''; CALL flush(nf)
    DO k=m1+1,2*m1; WRITE (nf,FMT='(E15.5,1X,$)') var3(k); END DO; WRITE (nf,FMT='(A)') ''; CALL flush(nf)
    WRITE (nf,FMT='(A)') ''; CALL flush(nf)
    WRITE (nf,FMT='(2(E25.15,1X))') cs1,cs2; CALL flush(nf)
   
  END SUBROUTINE dumpData

  SUBROUTINE packBuffer(m1,var1,var2,buffer)
     INTEGER,INTENT(IN)  :: m1
     REAL   ,INTENT(IN)  :: var1(m1)
     LOGICAL, INTENT(IN) :: var2(m1)
     REAL   ,INTENT(OUT) :: buffer(2*m1)
     
     INTEGER :: i
     REAL :: vaux
     REAL,PARAMETER :: r_true=7.77e+07
     
     DO i=1,m1
        buffer(i)=var1(i)
	IF(var2(i)) THEN
	   vaux=r_true
	ELSE
	   vaux=0.0
	END IF
	buffer(m1+i)=vaux
     END DO
     
  END SUBROUTINE packBuffer

  SUBROUTINE unPackBuffer(m1,buffer,var1,var2)
     INTEGER,INTENT(IN)  :: m1
     REAL   ,INTENT(IN)  :: buffer(2*m1)
     REAL   ,INTENT(OUT) :: var1(m1)
     LOGICAL,INTENT(OUT) :: var2(m1)
     
     INTEGER :: i
     
     
     DO i=1,m1
        var1(i)=buffer(i)
	IF(buffer(m1+i)>10.0) THEN
	   var2(i)=.TRUE.
	ELSE
	   var2(i)=.FALSE.
	END IF
     END DO
     
  END SUBROUTINE unPackBuffer



!-------------------------
SUBROUTINE Advec1d(idime,ix,q0,u,den0,den1,dt,dxx,dd0,qn,line)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: ix,line
  INTEGER,INTENT(IN)                   :: idime
  REAL   ,INTENT(in),DIMENSION(ix) :: q0
  REAL   ,INTENT(in),DIMENSION(ix)   :: u
  REAL   ,INTENT(in),DIMENSION(ix)     :: den0
  REAL   ,INTENT(in),DIMENSION(ix)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(ix)     :: dxx
  REAL   ,INTENT(in),DIMENSION(ix)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(ix)    :: qn

  REAL,DIMENSION(ix)                 :: flux
  REAL,DIMENSION(ix)                   :: vcmax
  REAL,DIMENSION(ix)                   :: vcmin
  LOGICAL,DIMENSION(ix)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1
     imxmn(i)=q0(i)>=max(q0(i-1),q0(i+1)) .or. & !=true if local
              q0(i)<=min(q0(i-1),q0(i+1))       !       extrema
     ck1= q0(i)
     ck2= q0(i)
     if(u(i  )< zr0) ck1= q0(i+1)
     if(u(i-1)>=zr0) ck2= q0(i-1)
     vcmax(i)= max( q0(i), ck1, ck2 )                      ! Eq-7
     vcmin(i)= min( q0(i), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1)>=zr0) flux(1)= q0(1)*u(1)*dt*dd0(1)
  DO i=2,idime-1
  
  !print*,'part 1'
     IF(u(i)<zr0) CYCLE
     IF(u(i-1)<zr0) THEN
        flux(i)= q0(i)*u(i)*dt*dd0(i)    !  outflow-only cell
     ELSE                              !      use upstream
        x1= dt*u(i)/dxx(i)               ! Courant number
        x1n= (1.-x1)*(q0(i+1)-q0(i-1))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(i) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(i-1)) cf= q0(i) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(i+1)) cf= q0(i) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(I) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(i),q0(i+1))  ), MAX(q0(i),q0(i+1)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(i)= MAX(  vcmin(i),   MIN(   vcmax(i),          &   !eq-3&8
        (q0(i)*den0(i)-x1*cf1*dd0(i)+flux(i-1)/dxx(i))/den1(i) ))

	! for debug purposes only
        if(debug==on)  QN(I)=(Q0(I)*DEN0(I)-X1*CF1*DD0(I)+FLUX(I-1)/DXX(I))/DEN1(I)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(i)= dxx(i)*(q0(i)*den0(i) - qn(i)*den1(i)) + flux(i-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  IF(u(idime-1)<zr0) flux(idime-1)=q0(idime)*u(idime-1)*dt*dd0(idime-1)
  DO i=idime-1,2,-1
 ! print*,'part 2'
       IF(u(i-1)>=zr0) THEN           ! Inflow-only cell
        IF(u(i)<zr0) qn(i)=  MAX(  vcmin(i),   MIN(   vcmax(i),&
           (q0(i)*den0(i)-flux(i)/dxx(i) + flux(i-1)/dxx(i))/den1(i) ))
        ELSE
           x1=  dt*ABS(u(i-1))/dxx(i)     ! Courant number
           x1n= (1.-x1)*(q0(i-1)-q0(i+1))/4.
           cf= q0(i) + x1n                                       !Eq-4b
           IF(imxmn(i+1)) cf= q0(i) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(i-1)) cf= q0(i) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(i),q0(i-1)) ), MAX(q0(i),q0(i-1)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(i)>=zr0) cf1= q0(i)     ! outflow-only cell upstream
           qn(i)= MAX(  vcmin(i),  MIN(   vcmax(i),           &   !Eq-3&8
               (q0(i)*den0(i)-flux(i)/dxx(i)-x1*cf1*dd0(i-1))/den1(i) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(I)=(Q0(I)*DEN0(I)-FLUX(I)/DXX(I)-X1*CF1*DD0(I-1))/DEN1(I)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(i-1)=dxx(i)*(qn(i)*den1(i) - q0(i)*den0(i)) + flux(i)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
  END SUBROUTINE Advec1d
!-------------------------------------------------------------------
 
  SUBROUTINE prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,u3d,v3d,w3d&
			    ,dxt,dyt &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )

   
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d     &
                                        ,dd0_3d ,dd0_3du &
					,dd0_3dv,dd0_3dw 
  
   !- local var
   REAL   :: dtlto2
   INTEGER :: i,j,k
   REAL  :: ai0s  =  25.0
   REAL  :: aj0s  =  50.0
   REAL  :: umx   =   10.0
   REAL,PARAMETER :: pii   =   3.141592653589793
   REAL    :: umax  =   0.0
   REAL    :: anrev,curnt,rx,xa,ilop,iwndty,nrec
   
   dtlto2 =  10.!*dtlt
    
     !WRITE(6,*) ' Wind fields?  0-rotating; or  1-divergent winds'
     iwndty = 0  ! 0-rotating
     !iwndty = 1  ! 1-divergent winds

     IF(iwndty==1) ai0s= 50.5
     ilop= ai0s-21.  ! needed for printouts
     ! Define wind fields (rotation or divergent) and initial mixing ratios
     !  Cone at (25,50) for rotating winds; Cone at (50,50) divergent winds
     DO k = 1,m1
   	DO  j=m3,1,-1
   	   DO  i=1,m2
   	   
   	       dd0_3d (k,i,j)=1.
	       dd0_3du(k,i,j)=1. 
               dd0_3dv(k,i,j)=1.
	       dd0_3dw(k,i,j)=1. 

   	       u3d(k,i,j)= -2.*umx*(REAL(j)-REAL(110)/2.-.5)/REAL(110)*dtlto2
   	       v3d(k,i,j)=  2.*umx*(REAL(i)-REAL(100)/2.-.5)/REAL(100)*dtlto2
 	       w3d(k,i,j)= 0.                                        *dtlto2
   	   
   	      IF(iwndty==1) THEN
   		 xa=pii/25.
   		 IF(J>0) u3d(k,i,j)=umx*SIN(xa*REAL(i))*SIN(xa*(REAL(j)))
   		 IF(I>0) v3d(k,i,j)=umx*COS(xa*(REAL(i)-.5))*COS(XA*(REAL(j)+.5))
   	      END IF
   	   
   	   
   	      umax= MAX(ABS(u3d(k,i,j)),ABS(v3d(k,i,j)),umax)
   	      rx= SQRT((REAL(i)-ai0s)**2.+(REAL(j)-aj0s)**2.)

   	   END DO !i
   	END DO !j
     END DO !k  

     !IF(iwndty==0) THEN
   	!WRITE(6,*) ' #iterations per revolution? must be>',pii*(m2-1)
   	!WRITE(6,*) '	Followed by the number of revolutions desired'
   	!srf- READ(5,*) NITPREV, ANREV
   	!nitprev=320; anrev=1
   	!dt= pii*REAL(m2-1)*dx/umax/REAL(nitprev)
   	!nstep= REAL(nitprev)*anrev + .9999

     !ELSE
   	!WRITE(6,*) ' Courant #?  eg. C=.6594 (to match Staniforth figs)'
   	!WRITE(6,*) '  Followed by # of iterations? e.g. 20,40,60,80,400'
   	! READ(5,*) CURNT, NSTEP
   	!curnt=0.6594 ;nstep=80
   	!dt= curnt*dx/umax
     !END IF

     !WRITE(6,*) ' Max wind speed, Courant number=',umax,umax*dt/dx
     !WRITE(6,*) nstep,' iterations required'
     !CALL FLUSH(6)


  END SUBROUTINE prepare_theor_winds
!LFR>  !----------------------------------------------------------------------------  
!LFR> 
!LFR>   SUBROUTINE check_mass_conservation
!LFR>   use mem_grid
!LFR>   use var_tables
!LFR>   use node_mod
!LFR>   use mem_basic
!LFR>   use mem_aer1, only : AEROSOL,num_scalar_aer_1st
!LFR>   use mem_chem1, only : NSPECIES_TRANSPORTED
!LFR>   
!LFR>   implicit none
!LFR>   real masscon, initalmass
!LFR>   integer :: mxyzp,n
!LFR>   mxyzp = mxp * myp * mzp
!LFR>   initalmass=6.8036989E+17
!LFR>   
!LFR>   !do n = 1,num_scalar(ngrid)
!LFR>   !srf - somente para gases, veja rotina original para todos os demais escalares
!LFR>   do n=num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED)+1,num_scalar(ngrid)
!LFR>      !------------------- TMP (SO para 1 primeiro gas)
!LFR>      if(n>num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED)+1) CYCLE ! TMP
!LFR>      !------------------- TMP
!LFR>   
!LFR>   
!LFR> !     call calc_masscon(mzp,mxp,myp,scalar_tab(n,ngrid)%var_p &
!LFR> !  		 ,grid_g(ngrid)%dxt    (1,1)	&
!LFR> !  		 ,grid_g(ngrid)%dyt    (1,1)	&
!LFR> !  		 ,basic_g(ngrid)%dn0 &
!LFR> !  		 , masscon)
!LFR>   
!LFR>   print*,'N masscon=',n, masscon,100.*(masscon-initalmass)/initalmass
!LFR> 
!LFR>   enddo
!LFR>   end SUBROUTINE check_mass_conservation
!LFR>   !-------------------------------------------------------

  SUBROUTINE calc_masscon(m1,m2,m3,var_p, dxt,dyt,dn0,masscon)
  use mem_grid, only : dzt
  implicit none
  integer, intent(in) :: m1,m2,m3
  real, dimension(m1,m2,m3), intent(in) :: var_p,dn0
  real, intent(out) :: masscon
  real, dimension(m2,m3), intent(in) :: dxt,dyt
  integer :: i,j,k
  real vol
  
  masscon=0.
  do j=2,m3-1; do i=2,m2-1
  
  
      do k=2,m1-1
    !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
    vol= 1./( dzt(k)*dxt(i,j)*dyt(i,j)) 
    masscon=masscon+var_p(k,i,j)*dn0(k,i,j)*vol

  enddo;enddo;enddo

  end SUBROUTINE calc_masscon
  
  
   SUBROUTINE InitCommMatrix(Mynum,master_num)
   
      use rpara, only: &
               nmachs,ixb,ixe,iyb,iye
      use mem_grid, only: &
               ngrids, nnxp,nnyp
!LFR>       use ParLib50, only: &
!LFR>                parf_bcast  !subroutine     
	       
     
     IMPLICIT NONE
      
      INCLUDE 'mpif.h'
  
      INTEGER, INTENT(IN) :: mynum,master_num
      
      INTEGER :: i,j,ngrid,nproc,nprocrec,fromto,ncount,nvar,ierr,inx,iny,ng
      INTEGER :: procmap1d(ngrids*maxval(nnxp)*maxval(nnyp))
      INTEGER :: info(1),ihandle(2),countItag
      !INTEGER :: inittag(nmachs)

      ALLOCATE(procMap(ngrids,maxval(nnxp),maxval(nnyp)));procMap=0
      ALLOCATE(handle(nmachs+100000))      
      
      IF(mynum==master_num) THEN !I am the master
        
         PRINT *,'Notice from the master:'
         PRINT *,'Number of processors: ',nmachs
         PRINT *,'Master node is the  : ',mynum
         PRINT *,'ngrids, maxval(nnxp),maxval(nnyp)=',ngrids, maxval(nnxp),maxval(nnyp)
         ALLOCATE(itag1(ngrids,nmachs,nmachs,nvariables))
         ALLOCATE(itag2(ngrids,nmachs,nmachs,nvariables))
 

	 info(1)=nmachs
         CALL flush(6)    

         DO ngrid=1,ngrids
            DO nproc=1,nmachs
               DO i=ixb(nproc,ngrid),ixe(nproc,ngrid)
                  DO j=iyb(nproc,ngrid),iye(nproc,ngrid)
	             procMap(ngrid,i,j)=nproc   
	          END DO
               END DO
	    END DO
         END DO
	
	  WRITE(*,FMT='(A)') '--(.)(.)-- Map of processors for advection --(.)(.)--';CALL flush(6)
	  DO i=1,nnxp(1)
	     WRITE(*,FMT='(A,I3.3,1X,$)') 'I= ',i
	     DO j=1,nnyp(1)
		WRITE(*,FMT='(I2.2,A,$)') procMap(1,i,j),' ';CALL flush(6)
	     END DO
	     WRITE(*,FMT='(A)') '';CALL flush(6)
	  END DO

         ncount=0
	 DO ngrid=1,ngrids
	    DO i=1,maxval(nnxp)
	       DO j=1,maxval(nnyp)
	             ncount=ncount+1
	             procmap1d(ncount)=procMap(ngrid,i,j)
	       END DO
	    END DO
	 END DO
	 DO nproc=1,nmachs	
	    !PRINT *,'Notice from ',mynum,' procmap sended to proc #',nproc; CALL flush(6)    
	    CALL send(info,1,nproc,nproc+50000)
	    CALL send(procmap1d,ngrids*maxval(nnxp)*maxval(nnyp),nproc,nproc)
	 END DO
      ELSE !I am a slave      
	 CALL get(info,1,0,mynum+50000)
	 CALL get(procmap1d,ngrids*maxval(nnxp)*maxval(nnyp),master_num,mynum)
	 CALL parf_wait_nostatus(handle(mynum+50000))
	 CALL parf_wait_nostatus(handle(mynum))
         !WRITE (mynum+90,FMT='(A)') 'Notice from slave   :'
         !WRITE (mynum+90,FMT='(A,I3.3)') 'Number of processors: ',info(1)
         !WRITE (mynum+90,FMT='(A,I3.3)') 'Master node is the  : ',mynum
         !WRITE (mynum+90,FMT='(A,3(I3.3,1X))') 'ngrids, maxval(nnxp),maxval(nnyp)=',ngrids, maxval(nnxp),maxval(nnyp)
         !CALL flush(90+mynum)
         ALLOCATE(itag1(ngrids,info(1),info(1),nvariables))
         ALLOCATE(itag2(ngrids,info(1),info(1),nvariables))
         numberOfProcessors=info(1)
         ncount=0
         DO ngrid=1,ngrids
	    DO i=1,maxval(nnxp)
	       DO j=1,maxval(nnyp)
	             ncount=ncount+1
	             procMap(ngrid,i,j)=procmap1d(ncount)
	          END DO
	    END DO
         END DO      
      END IF

      ncount=1
      DO ngrid=1,ngrids
         DO nproc=1,info(1)
             DO nprocrec=1,info(1)
	        DO nvar=1,nvariables
		         itag1(ngrid,nproc,nprocrec,nvar)=ncount
		         itag2(ngrid,nproc,nprocrec,nvar)=ncount+1
		END DO
	     END DO
	  END DO
      END DO 
      !PRINT *,'ncount=',ncount;CALL flush(6)
	 
      ALLOCATE (west(ngrids, info(1),maxval(nnyp)),east(ngrids,info(1),maxval(nnxp)), &
		north(ngrids,info(1),maxval(nnxp)),south(ngrids,info(1),maxval(nnyp)))
      ALLOCATE(lineTag(ngrids,info(1),info(1),maxval(nnxp),2))
      ALLOCATE(columnTag(ngrids,info(1),info(1),maxval(nnyp),2))
      lineTag=0
      columnTag=0

      CountItag=0
      west=0;east=0;north=0;south=0
      DO ng=1,ngrids
         DO i=2,nnxp(ng)-1
            DO j=2,nnyp(ng)-1
	       countItag=countItag+1
               IF (procMap(ng,i,j)/=procMap(ng,i,j-1) .and. procMap(ng,i,j-1)/=0) THEN
	          west(ng,procMap(ng,i,j),i)=procMap(ng,i,j-1)
	          east(ng,procMap(ng,i,j-1),i)=procMap(ng,i,j)
	          LineTag(ng,procMap(ng,i,j-1),procMap(ng,i,j),i,1)=countItag
	          LineTag(ng,procMap(ng,i,j),procMap(ng,i,j-1),i,2)=countItag	   
	          countItag=countItag+1		  
	          LineTag(ng,procMap(ng,i,j-1),procMap(ng,i,j),i,2)=countItag
	          LineTag(ng,procMap(ng,i,j),procMap(ng,i,j-1),i,1)=countItag	                  
	       END IF
            END DO
         END DO
         DO j=2,nnyp(ng)-1
            DO i=2,nnxp(ng)-1
               countItag=countItag+1
               IF (procMap(ng,i,j)/=procMap(ng,i-1,j) .and. procMap(ng,i-1,j)/=0) THEN
	          north(ng,procMap(ng,i,j),j)=procMap(ng,i-1,j)
	          south(ng,procMap(ng,i-1,j),j)=procMap(ng,i,j)
	          columnTag(ng,procMap(ng,i,j),procMap(ng,i-1,j),j,1)=countItag
	          columnTag(ng,procMap(ng,i-1,j),procMap(ng,i,j),j,2)=countItag
	          countItag=countItag+1		  
	          columnTag(ng,procMap(ng,i,j),procMap(ng,i-1,j),j,2)=countItag
	          columnTag(ng,procMap(ng,i-1,j),procMap(ng,i,j),j,1)=countItag
	      END IF
            END DO
         END DO
      END DO
      
      DEALLOCATE(handle)
      ALLOCATE(handle(countItag))

      ALLOCATE(comm(ngrids,numberOfProcessors,4,2,numberOfProcessors)) !Allocating the total of processors and 4 directions (N,S,E,W)

      

      IF(mynum==2) THEN
      DO nproc=1,info(1)
	DO ng=1,info(1)
	   DO i=2,nnxp(1)
	      IF(lineTag(1,nproc,ng,i,1)/=0 .or. lineTag(1,nproc,ng,i,2)/=0) &
	      WRITE (*,FMT='(5(I5.5,2X))') nproc,ng,i,lineTag(1,nproc,ng,i,1),lineTag(1,nproc,ng,i,2)
	      IF(columnTag(1,nproc,ng,i,1)/=0 .or.  columnTag(1,nproc,ng,i,2)/=0) &
	      WRITE (*,FMT='(5(I5.5,2X))') nproc,ng,i,ColumnTag(1,nproc,ng,i,1),ColumnTag(1,nproc,ng,i,2)
	   END DO
	END DO
      END DO
      END IF	
     !STOP 500    
      
     !Just for test
     OPEN(27,file='walcekInit',form='unformatted',STATUS='replace')
     WRITE(27) ngrids, info(1),maxval(nnyp),maxval(nnxp),countItag
     WRITE(27) north,south,east,west
     WRITE(27) LineTag,columnTag
     CLOSE(27)
      
     call create_watch(w, name=(/ 'advct X', 'waiting', 'big lo ','advec T'/)) 
     call create_watch(w1, name=(/'Wait Up', 'Wait Dn', 'process','Sending', 'initial','fase 6'/))     
     iounit=80+mynum
     
   END SUBROUTINE InitCommMatrix  
  !-------------------------------------------------------
   SUBROUTINE prepareCommunication(ia,iz,ja,jz)
      INTEGER,INTENT(IN) :: ia,iz,ja,jz
      INTEGER,PARAMETER :: iNorth=1,iSouth=2,iEast=3,iWest=4
      INTEGER,ALLOCATABLE,DIMENSION(:) :: iCount,firstPoint,lastPoint
      
      INTEGER :: i,j,ng,ii,jj

 !    DO ng=1,ngrids     !grids
 !       DO np=1,info(1) !num. of  proc.
 !	   DO dir=1,4   !N,S,E,W
 !	      DO s=1,2  !from,to
 
      !PRINT *,mynum,ia,iz,ja,jz,nodei0(mynum,ngrid),nodej0(mynum,ngrid),numberOfProcessors
      ALLOCATE(iCount(numberOfProcessors))
      ALLOCATE(firstPoint(numberOfProcessors))
      ALLOCATE(lastPoint(numberOfProcessors))
      icount=0
      firstPoint=iz
      lastPoint=ia
      !West
      j=ja+nodej0(mynum,ngrid)
      DO ng=1,ngrid
         DO i=ia,iz
	    ii=i+nodei0(mynum,ngrid)
	    IF(procMap(ng,ii,j)/=procMap(ng,ii,j-1) .and. procMap(ng,ii,j-1)/=0) THEN
	       icount(procMap(ng,ii,j-1))=icount(procMap(ng,ii,j-1))+1
	       firstPoint(procMap(ng,ii,j-1))=min(firstPoint(procMap(ng,ii,j-1)),i)
	       lastPoint(procMap(ng,ii,j-1))=max(lastPoint(procMap(ng,ii,j-1)),i)	    
	    END IF
	 END DO
	 DO i=1,numberOfProcessors
            IF(icount(i)/=0) THEN
               comm(ng,mynum,iWest,fromP,i)%npoints=icount(i)
	       comm(ng,mynum,iWest,fromP,i)%iniP=firstPoint(i)
	       comm(ng,mynum,iWest,fromP,i)%endP=lastPoint(i)
	       PRINT *,"West: ",mynum,i,&
	         comm(ng,mynum,iWest,fromP,i)%npoints,comm(ng,mynum,iWest,fromP,i)%iniP, &
		 comm(ng,mynum,iWest,fromP,i)%endP;CALL flush(6)
	    END IF
         END DO
      END DO
      icount=0
      firstPoint=iz
      lastPoint=ia
      !east
      j=jz+nodej0(mynum,ngrid)
      DO ng=1,ngrid
         DO i=ia,iz
	    ii=i+nodei0(mynum,ngrid)
	    IF(procMap(ng,ii,j)/=procMap(ng,ii,j+1) .and. procMap(ng,ii,j+1)/=0) THEN
	       icount(procMap(ng,ii,j+1))=icount(procMap(ng,ii,j+1))+1
	       firstPoint(procMap(ng,ii,j+1))=min(firstPoint(procMap(ng,ii,j+1)),i)
	       lastPoint(procMap(ng,ii,j+1))=max(lastPoint(procMap(ng,ii,j+1)),i)	    
	    END IF
	 END DO
	 DO i=1,numberOfProcessors
            IF(icount(i)/=0) THEN
               comm(ng,mynum,ieast,fromP,i)%npoints=icount(i)
	       comm(ng,mynum,ieast,fromP,i)%iniP=firstPoint(i)
	       comm(ng,mynum,ieast,fromP,i)%endP=lastPoint(i)
	       PRINT *,"East: ",mynum,i,&
	         comm(ng,mynum,ieast,fromP,i)%npoints,comm(ng,mynum,ieast,fromP,i)%iniP, &
		 comm(ng,mynum,ieast,fromP,i)%endP;CALL flush(6)
	    END IF
         END DO
      END DO
      icount=0
      firstPoint=jz
      lastPoint=ja
      !south
      i=iz+nodei0(mynum,ngrid)
      DO ng=1,ngrid
         DO j=ja,jz
	    jj=j+nodej0(mynum,ngrid)
	    IF(procMap(ng,i,jj)/=procMap(ng,i+1,jj) .and. procMap(ng,i+1,jj)/=0) THEN
	       icount(procMap(ng,i+1,jj))=icount(procMap(ng,i+1,jj))+1
	       firstPoint(procMap(ng,i+1,jj))=min(firstPoint(procMap(ng,i+1,jj)),j)
	       lastPoint(procMap(ng,i+1,jj))=max(lastPoint(procMap(ng,i+1,jj)),j)	    
	    END IF
	 END DO
	 DO i=1,numberOfProcessors
            IF(icount(i)/=0) THEN
               comm(ng,mynum,isouth,fromP,i)%npoints=icount(i)
	       comm(ng,mynum,isouth,fromP,i)%iniP=firstPoint(i)
	       comm(ng,mynum,isouth,fromP,i)%endP=lastPoint(i)
	       PRINT *,"South ",mynum,i,&
	         comm(ng,mynum,isouth,fromP,i)%npoints,comm(ng,mynum,isouth,fromP,i)%iniP, &
		 comm(ng,mynum,isouth,fromP,i)%endP;CALL flush(6)
	    END IF
         END DO
      END DO
       icount=0
      firstPoint=jz
      lastPoint=ja
      !North
      i=ia+nodei0(mynum,ngrid)
      DO ng=1,ngrid
         DO j=ja,jz
	    jj=j+nodej0(mynum,ngrid)
	    IF(procMap(ng,i,jj)/=procMap(ng,i-1,jj) .and. procMap(ng,i-1,jj)/=0) THEN
	       icount(procMap(ng,i-1,jj))=icount(procMap(ng,i-1,jj))+1
	       firstPoint(procMap(ng,i-1,jj))=min(firstPoint(procMap(ng,i-1,jj)),j)
	       lastPoint(procMap(ng,i-1,jj))=max(lastPoint(procMap(ng,i-1,jj)),j)	    
	    END IF
	 END DO
	 DO i=1,numberOfProcessors
            IF(icount(i)/=0) THEN
               comm(ng,mynum,isouth,fromP,i)%npoints=icount(i)
	       comm(ng,mynum,isouth,fromP,i)%iniP=firstPoint(i)
	       comm(ng,mynum,isouth,fromP,i)%endP=lastPoint(i)
	       PRINT *,"North ",mynum,i,&
	         comm(ng,mynum,isouth,fromP,i)%npoints,comm(ng,mynum,isouth,fromP,i)%iniP, &
		 comm(ng,mynum,isouth,fromP,i)%endP;CALL flush(6)
	    END IF
         END DO
      END DO
     


   END SUBROUTINE prepareCommunication
  !-------------------------------------------------------
  
  
  SUBROUTINE get_real(var,m1,from,tag)
       
     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     REAL,INTENT(OUT)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_real(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_real
  
  SUBROUTINE get_int(var,m1,from,tag)
       
     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     INTEGER,INTENT(OUT)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_int(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_int

  SUBROUTINE get_logical(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     LOGICAL,INTENT(OUT) :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_logical(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_logical
  
  SUBROUTINE send_real(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     REAL,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_real(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_real  
  
  SUBROUTINE send_int(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     INTEGER,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_int(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_int 

  SUBROUTINE send_logical(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     LOGICAL,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_logical(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_logical 
  

!---------------------------------------------------------------------------

!LFR> !----------------------------------------------------
!LFR>  SUBROUTINE advect_mnt_stalone()
!LFR>  
!LFR>    use node_mod, only : nodei0,nodej0
!LFR>    
!LFR>   IMPLICIT NONE
!LFR>   INTEGER :: m1
!LFR>   INTEGER :: m2
!LFR>   INTEGER :: m3
!LFR>   INTEGER :: ia
!LFR>   INTEGER :: iz
!LFR>   INTEGER :: ja
!LFR>   INTEGER :: jz,n
!LFR>   INTEGER :: mynum
!LFR>   
!LFR>   REAL :: dt
!LFR>   
!LFR>   REAL,DIMENSION(:)     , ALLOCATABLE  :: dzt
!LFR>   REAL,DIMENSION(:,:)   , ALLOCATABLE  :: dxt,dyt,rtgt
!LFR>   
!LFR>   REAL,DIMENSION(:),       ALLOCATABLE :: dztW 
!LFR>   REAL,DIMENSION(:,:)   , ALLOCATABLE :: dxtW,dytW
!LFR> 
!LFR>   REAL,DIMENSION(:,:,:), ALLOCATABLE :: vc3d_in! initial mixing ratio 
!LFR>   REAL,DIMENSION(:,:,:), ALLOCATABLE :: u3d,v3d,w3d
!LFR>   REAL,DIMENSION(:,:,:), ALLOCATABLE :: dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d
!LFR> 
!LFR>   REAL,DIMENSION(:,:,:),ALLOCATABLE :: vc3d_x,vc3d_y
!LFR>   REAL,DIMENSION(:,:,:), ALLOCATABLE :: vc3d_out ! updated mixing ratio due advection
!LFR> 
!LFR>   !- local var
!LFR>   REAL,DIMENSION(:),ALLOCATABLE               :: dxx  
!LFR>   REAL,DIMENSION(:,:),ALLOCATABLE            :: dxy  
!LFR>   real masscon,initialmass,vol
!LFR>   integer :: i,j,k,nrec,it,nmachs,ngrids
!LFR>   INTEGER :: step
!LFR>   
!LFR>   CHARACTER(LEN=2) :: nf
!LFR> 
!LFR>   DO step=1,20
!LFR> 
!LFR>   !For test standalone write
!LFR>   WRITE(nf,FMT='(I2.2)') step
!LFR>   !OPEN(27,file='walcekDataIn.'//nf,form='unformatted',STATUS='replace')
!LFR>   !WRITE(27) m1,m2,m3,ia,iz,ja,jz
!LFR>   !WRITE(27) dzt,dxt,dyt,rtgt,dztW,dxtW,dytW
!LFR>   !WRITE(27) vc3d_in,u3d,v3d,w3d,dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d
!LFR>   !PRINT *,'Done!'
!LFR>   !CLOSE(27)
!LFR>   OPEN(27,file='walcekDataIn.'//nf,form='unformatted',STATUS='old')
!LFR>   READ(27) m1,m2,m3,ia,iz,ja,jz
!LFR>   IF(step==1) THEN
!LFR>      ALLOCATE(dzt(m1),dztW(m1),dxx(m1))
!LFR>      ALLOCATE(dxt(m2,m3),dyt(m2,m3),rtgt(m2,m3),dxtW(m2,m3),dytW(m2,m3),dxy(m2,m3))
!LFR>      ALLOCATE(vc3d_in(m1,m2,m3),u3d(m1,m2,m3),v3d(m1,m2,m3),w3d(m1,m2,m3))
!LFR>      ALLOCATE(dd0_3d(m1,m2,m3),den0_3d(m1,m2,m3),den1_3d(m1,m2,m3),den2_3d(m1,m2,m3))
!LFR>      ALLOCATE(den3_3d(m1,m2,m3),vc3d_x(m1,m2,m3),vc3d_y(m1,m2,m3),vc3d_out(m1,m2,m3))
!LFR>   END IF
!LFR>   READ(27) dt,nodei0,nodej0
!LFR>   READ(27) dzt,dxt,dyt,rtgt,dztW,dxtW,dytW
!LFR>   READ(27) vc3d_in,u3d,v3d,w3d,dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d
!LFR>   PRINT *,'Done!'
!LFR>   CLOSE(27)
!LFR>   !End Test read state
!LFR>   
!LFR>   ! Do I-advection first over all J-k lines
!LFR>  
!LFR>   call Advec3d_X(m1,m2,m3,ia,iz,ja,jz,vc3d_in,u3d,den0_3d, &
!LFR>                  den1_3d,dt,dxtW,dd0_3d,vc3d_x,mynum,n )
!LFR> 
!LFR>   
!LFR>   !- do J-advection next over all I-k lines
!LFR> 
!LFR>   call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3d_x ,v3d,den1_3d, &
!LFR>                  den2_3d,dt,dytW,dd0_3d,vc3d_y )
!LFR>  
!LFR>   
!LFR>   !- do k-advection next over all I-J lines
!LFR>   ! vc3dout_Z = razao de mistura final (apos a adveccao nas 3 dimensoes)
!LFR>   call Advec3d_Z(m1,m2,m3,ia,iz,ja,jz,vc3d_y,w3d,den2_3d, &
!LFR>                  den3_3d,dt,dztW,dd0_3d,vc3d_out)
!LFR>   
!LFR>   
!LFR>    !For test standalone
!LFR>   OPEN(19,file='walcekDataOut.'//nf//'.gra',         &
!LFR>        form='unformatted',access='direct',status='unknown',  &
!LFR>        recl=4*m2*m3)
!LFR>   nrec =0
!LFR>   do k=1,m1
!LFR>     nrec=nrec+1
!LFR>     write (19,rec=nrec) vc3d_out(k,:,:)
!LFR>   enddo
!LFR>   CLOSE(19)
!LFR>  !End test
!LFR> 
!LFR>   END DO
!LFR>   STOP 'end of test'
!LFR>   !-------------------------------
!LFR>   !- check mass conservation (after advection)
!LFR>   !masscon=0.
!LFR>   !do j=2,m3-1; do i=2,m2-1
!LFR>   !  do k=2,m1-1
!LFR>   !  !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
!LFR>   !  vol= 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j) ! true volume
!LFR>   !  masscon=masscon+vc3dout_Z(k,i,j)*dd0_3d(k,i,j)*vol
!LFR>   !
!LFR>   !enddo;enddo;enddo
!LFR>   !print*,'masscon=',initialmass,masscon,100.*(masscon-initialmass)/initialmass
!LFR>   !print*,'MAXVALS=', maxval(vc3din)*28.96/28, maxval(vc3dout_Z)*28.96/28
!LFR>   !call flush(6)
!LFR>   !-------------------------------
!LFR> 
!LFR> !---------------------------------------------------------------------
!LFR> !- for theoretical experiments
!LFR> !---------------------------------------------------------------------
!LFR> !  vc3din=vc3dout_Z
!LFR> !
!LFR> ! ENDDO
!LFR> ! close (19)
!LFR> ! stop 333
!LFR> !---------------------------------------------------------------------
!LFR> !---------------------------------------------------------------------
!LFR>    
!LFR>  END SUBROUTINE advect_mnt_stalone



END MODULE monotonic_adv
